# Appendix Figures
source("analysis.r")

# -----------------------------------
# cov balance plot -- Figure A.1
# -----------------------------------

# recode variable names
recode_var <- c(
  "years_emp_1" = "Years employed: < 5 yrs",
  "position_1" = "Top - Business Headquarters Class",
  "position_2" = "Position: Manager - Chief class",
  "position_3" = "Position: Chief class and below",
  "college_1" = "College degree",
  "manufacturing_1" = "Manufacturing",
  "employee_2" = "Employee: 300-5000",
  "employee_3" = "Employee: > 5000",
  "TKY_OSK_1" = "Location: Tokyo and Osaka",
  "estab_year_1" = "Established year: pre-1945",
  "foreign_owned_1" = "Japanese ownership",
  "sales_2" = "Sales: > 50B and < 100B yen",
  "sales_3" = "Sales: > 100B yen",
  "capital_2" = "Capital: > 1M and < 10B yen",
  "capital_3" = "Capital: > 10B yen",
  "sanction_impact" = "Sanction impact: negative",
  "sanction_second" = "Support secondary sanction",
  "sanction" = "Support sanction", 
  "income_1" = "Income: 0-6M yen",
  "age_2" = "Age: 45-60 years old",
  "age_3" = "Age: 60+ years old"
)


order_var = c(  "employee_2",
                "employee_3",
                "TKY_OSK_1",
                "manufacturing_1",
                "estab_year_1",
                "foreign_owned_1",
                "capital_2",
                "capital_3",
                "sales_2",
                "sales_3",
                "position_1",
                "position_2",
                "position_3",
                "years_emp_1",
                "age_2",
                "age_3",
                "income_1",
                "college_1",
                "sanction_impact",
                "sanction_second",
                "sanction")          

# covariates
covariate_ind = c("years_emp", "position", "college", "age", "income")
covariate_firm = c("manufacturing", "employee", "TKY_OSK", "estab_year", "foreign_owned", "capital", "sales")
pretreat_covs = u_pretreat_covs

u_sample %>% 
  dplyr::select(starts_with("u_treat")) %>% names -> treat_vec 

u_sample %>% 
  dplyr::select(
    all_of(c(
      "sanction","sanction_impact", "sanction_second",
      covariate_ind, covariate_firm, treat_vec))) %>% 
  fastDummies::dummy_cols(
    ., remove_selected_columns = TRUE,
    remove_first_dummy = TRUE) %>%
  dplyr::select(-contains("_NA")) %>%
  mutate(treatment = 1-u_treatment_control) -> balance_data


# compute the covariate balance 
var_remove <- c(treat_vec, "treatment")

ib_unmatch <- covariate_balance_atc(
  x_treat   = balance_data %>% 
    filter(treatment == 1) %>% 
    na.omit() %>%
    dplyr::select(-one_of(var_remove)),
  x_control = balance_data %>% 
    filter(treatment == 0) %>% 
    na.omit() %>%
    dplyr::select(-one_of(var_remove))
)

# plot
ib_unmatch %>%
  mutate(
    Variable = factor(Variable, levels = rev(order_var)),
    Variable = recode(Variable, !!!recode_var)) %>%
  filter(Variable != "NA") %>% 
  ggplot(., aes(x = Variable, y = `Standardized Difference`)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_point(size = 3) +
  coord_flip() +
  my_theme() +
  ylim(-1, 1) +
  xlab("") +
  ylab("Standardized Mean Difference") +
  geom_hline(yintercept = 0.25, color = '#6f6f6f', linetype = "dashed") +
  geom_hline(yintercept = -0.25, color = '#6f6f6f', linetype = "dashed")    

ggsave("plot/ukraine_cov_balance.pdf", width = 12, height = 14)



# ----------------------
# Russia business -- Figure D.1
# ----------------------

covariate_ind = c("position", "income")
covariate_firm = c("TKY_OSK", "capital")
pretreat_covs = c("sanction_impact")

df_russia <- u_sample %>% 
  filter(trade_russia_binary == 1) %>%
  mutate(u_stop = as.numeric(u_stop))

collect_models_w = c()

for(i in 1:3){
  treatment = treatment_vec[i]
  control = control_vec[i]

  run_lm(
    data = df_russia,
    treatment = treatment,
    control = control,
    outcome = "u_stop",
    covs = c(pretreat_covs, covariate_firm, covariate_ind),
    covs_include = TRUE) -> model_w
      
  collect_models_w[[i]] = model_w
}

# get results
data_fig4 <- bind_rows(
  collect_models_w[[3]] %>%
    get_data("US withdrawal"),
  collect_models_w[[2]] %>%
    get_data("Multiple withdrawal"),
  collect_models_w[[1]] %>%
    get_data("China stays")
)

data_fig4 %>%
  mutate(
      group  = c(rep(
          "Peer Conformity", 2), "Market Competition") %>% as_factor,
      label = factor(label, levels = c( "China stays", "US withdrawal", "Multiple withdrawal")),
      conf.low.95 = estimate - std.error * 1.96,
      conf.high.95 = estimate + std.error * 1.96,
      .width = 0.95,
      conf.low.90 = estimate - std.error * 1.64,
      conf.high.90 = estimate + std.error * 1.64,
      .width = 0.9) %>%
plot_ols_estimates(., "label") +
my_theme() +
ggtitle("Subset to Firms with Business in Russia:\nEffects on Support for Withdrawal from Russia") +
scale_x_continuous(breaks = c(-0.5, 0, 0.5)) +
theme(legend.position = "right") +
xlab("Changes in Support for Withdrawal from Russia") 


ggsave("plot/ukraine_H1-2_russia_ols.pdf", width = 8, height = 8/1.618)

# -----------------------------------
# impact of Ukraine War -- Figure D.2
# -----------------------------------

covariate_ind = c("years_emp", "position", "college")
covariate_firm = c("TKY_OSK", "estab_year", "foreign_owned")
pretreat_covs = c("sanction", "sanction_second")

treatment_vec = c("u_treatment_firm_china", "u_treatment_firm_multiple", "u_treatment_firm_us")

control_vec = c("u_treatment_firm_multiple", "u_treatment_control", "u_treatment_control")

# Positively affected or Not affected
df_impact <- u_sample %>% filter(sanction_impact == 0)

collect_models_impact = c()

for(i in 1:3){
  treatment = treatment_vec[i]
  control = control_vec[i]

 run_lm(
    data = df_impact,
    treatment = treatment,
    control = control,
    outcome = "u_stop",
    covs = c(pretreat_covs, covariate_firm, covariate_ind),
    covs_include = TRUE) -> model_impact

  collect_models_impact[[i]] = model_impact

}


# get results
data_impact <- bind_rows(
  collect_models_impact[[3]] %>%
    get_data("US withdrawal"),
  collect_models_impact[[2]] %>%
    get_data("Multiple withdrawal"),
  collect_models_impact[[1]] %>%
    get_data("China stays")) %>%
    mutate(new_label = "Positively affected\n or Not affected")


# negatively affected
df_impact_neg <- u_sample %>% filter(sanction_impact == 1)

collect_models_impact_neg = c()

for(i in 1:3){
  treatment = treatment_vec[i]
  control = control_vec[i]

  run_lm(
    data = df_impact_neg,
    treatment = treatment,
    control = control,
    outcome = "u_stop",
    covs = c(pretreat_covs, covariate_firm, covariate_ind),
    covs_include = TRUE) -> model_impact    

  collect_models_impact_neg[[i]] = model_impact

}

# get results
data_impact_neg <- bind_rows(
  collect_models_impact_neg[[3]] %>%
    get_data("US withdrawal"),
  collect_models_impact_neg[[2]] %>%
    get_data("Multiple withdrawal"),
  collect_models_impact_neg[[1]] %>%
    get_data("China stays")) %>%
    mutate(new_label = "Negatively affected")


# plot
bind_rows(data_impact, data_impact_neg) %>%
  mutate(
      label = factor(label, levels = c("China stays", "US withdrawal", "Multiple withdrawal")),
      group  = as_factor(new_label),
      conf.low.95 = estimate - std.error * 1.96,
      conf.high.95 = estimate + std.error * 1.96,
      .width = 0.95,
      conf.low.90 = estimate - std.error * 1.64,
      conf.high.90 = estimate + std.error * 1.64,
      .width = 0.9) %>%
    plot_ols_estimates(., "label") +
    my_theme() +
   ggtitle("")  +
    theme(legend.position = "right") +
    guides(color = guide_legend(reverse = TRUE)) 


ggsave("plot/ukraine_impact.pdf", width = 8, height = 6)


# ----------------------
# firm size -- Figure D.3
# ----------------------

firm_size <- unique(u_sample$employee) %>% na.omit()
firm_size <- as.vector(firm_size)

covariate_ind = c("years_emp", "college", "income")
covariate_firm = c("TKY_OSK", "estab_year", "foreign_owned")
pretreat_covs = c("sanction_impact")

firmsize_models_w = c()

for (k in c(1:length(firm_size))){
  
  df_industry <- u_sample %>% filter(employee == firm_size[[k]])

  collect_models_size = c()
    
  for(i in 1:3){
    treatment = treatment_vec[i]
    control = control_vec[i]
    
    run_lm(
        data = df_industry,
        treatment = treatment,
        control = control,
        outcome = "u_stop",
        covs = c(pretreat_covs, covariate_firm, covariate_ind),
        covs_include = TRUE) -> model_size
    
    collect_models_size[[i]] = model_size
    
  }

  # get results
  bind_rows(
    collect_models_size[[3]] %>%
        get_data("US withdrawal"),
    collect_models_size[[2]] %>%
        get_data("Multiple withdrawal"),
    collect_models_size[[1]] %>%
        get_data("China stays")) %>%
    mutate(firmsize = firm_size[k]) -> collect_size
    
  firmsize_models_w[[k]] = collect_size
  
}


est_firmsize_comb <- do.call(rbind, firmsize_models_w)

est_firmsize_comb$firmsize <- factor(est_firmsize_comb$firmsize,levels = c("1", "2", "3"),label = c( "Employee < 300","Employee between\n300 and 5000", "Employee > 5000"))

est_firmsize_comb %>% 
  mutate(
      label = factor(label, levels = c("China stays", "US withdrawal", "Multiple withdrawal")),
      group  = as_factor(firmsize),
      conf.low.95 = estimate - std.error * 1.96,
      conf.high.95 = estimate + std.error * 1.96,
      .width = 0.95,
      conf.low.90 = estimate - std.error * 1.64,
      conf.high.90 = estimate + std.error * 1.64,
      .width = 0.9) %>%
    plot_ols_estimates(., "label") +
    my_theme() +
  xlab("Estimates") + 
  ggtitle("") +
  scale_color_manual(
    values = c("#C93312","#aaaaaa","#006ad5"),
    guide = guide_legend(reverse = TRUE,title = "Employee")) +
  geom_vline(xintercept = 0, color = '#6f6f6f', linetype = "dashed") +
  ylab("") + 
  theme(legend.position = "right") 


ggsave("plot/ukraine_firmsize.pdf", width = 8, height = 6)



# ----------------------
# firm foreign business activities -- Figure D.4
# ----------------------

u_sample %>% 
  mutate(indirect_export_bin = ifelse(indirect_export == 1, 1, 0),
         indirect_import_bin = ifelse(indirect_import == 1, 1, 0),
         importer = ifelse(str_detect(trade_none, "1"), 0, 1) %>% replace_na(.,0),
         exporter = ifelse(str_detect(trade_none, "2"), 0, 1) %>% replace_na(.,0),
         foreign_subsidiary = ifelse(str_detect(trade_none, "3"), 0, 1) %>% replace_na(.,0),
         outsource = ifelse(str_detect(trade_none, "4"), 0, 1) %>% replace_na(.,0)) -> u_sample

covariate_ind = c("position", "college", "income")
covariate_firm = c("employee", "TKY_OSK")
pretreat_covs = c("sanction_impact")

product_vec <- c(
  "indirect_export_bin", "indirect_import_bin",
  "importer", "exporter", "foreign_subsidiary",
  "outsource")

firmtype_models_w = c()

for (k in c(1:length(product_vec))){
  
  df_firmtype <- u_sample %>% filter(!!sym(product_vec[k]) == 1)
  
  collect_models_bus = c()
    
  for(i in 1:3){
    treatment = treatment_vec[i]
    control = control_vec[i]
    
    run_lm(
        data = df_firmtype,
        treatment = treatment,
        control = control,
        outcome = "u_stop",
        covs = c(pretreat_covs, covariate_firm, covariate_ind),
        covs_include = TRUE) -> model_bus
    
    collect_models_bus[[i]] = model_bus
    
  }

  # get results
  bind_rows(
    collect_models_bus[[3]] %>%
        get_data("US withdrawal"),
    collect_models_bus[[2]] %>%
        get_data("Multiple withdrawal"),
    collect_models_bus[[1]] %>%
        get_data("China stays")) %>%
    mutate(firmtype = product_vec[k]) -> collect_bus
    
  firmtype_models_w[[k]] = collect_bus  
  
}


est_firmtype_comb <- do.call(rbind, firmtype_models_w)

est_firmtype_comb$firmtype <- factor(
    est_firmtype_comb$firmtype,levels = c("importer","exporter","indirect_import_bin",  "indirect_export_bin","foreign_subsidiary","outsource"),
    label = c(  "Importer", "Exporter","Indirect importer",  "Indirect exporter","Foreign subsidiary","Outsource"))

est_firmtype_comb %>% 
  mutate(
      label = factor(label, levels = c("China stays", "US withdrawal", "Multiple withdrawal")),
      group  = as_factor(firmtype),
      conf.low.95 = estimate - std.error * 1.96,
      conf.high.95 = estimate + std.error * 1.96,
      .width = 0.95,
      conf.low.90 = estimate - std.error * 1.64,
      conf.high.90 = estimate + std.error * 1.64,
      .width = 0.9) %>%
    plot_ols_estimates(., "label") +
    my_theme() +
  xlab("Estimates") + 
  ggtitle("") +
  scale_color_manual(values = c("#C93312","#aaaaaa","#006ad5", "#E69F00", "#56B4E9", "#009E73","#F0E442", "#0072B2", "#D55E00", "#CC79A7"),
                     guide = guide_legend(reverse = TRUE,
                                          title = "Employee")) +
  geom_vline(xintercept = 0, color = '#6f6f6f', linetype = "dashed") +
  ylab("") + 
  theme(legend.position = "right") 


ggsave("plot/ukraine_firmtype.pdf", width = 8, height = 8)



# -----------------------------------
# higher manager sample -- Figure D.5
# -----------------------------------

covariate_ind = c("years_emp", "college", "income")
covariate_firm = c("employee", "TKY_OSK", "foreign_owned")
pretreat_covs = c("sanction")

df_manager <- u_sample %>% filter(position == 1)
collect_models_ceo = c()

for(i in 1:3){
  treatment = treatment_vec[i]
  control = control_vec[i]

  run_lm(
    data = df_manager,
    treatment = treatment,
    control = control,
    outcome = "u_stop",
    covs = c(pretreat_covs, covariate_firm, covariate_ind),
    covs_include = TRUE) -> model_ceo    

  collect_models_ceo[[i]] = model_ceo

}

# get results
bind_rows(
    collect_models_ceo[[3]] %>%
        get_data("US withdrawal"),
    collect_models_ceo[[2]] %>%
        get_data("Multiple withdrawal"),
    collect_models_ceo[[1]] %>%
        get_data("China stays")) %>%
mutate(
    label = factor(label, levels = c("China stays", "US withdrawal", "Multiple withdrawal")),
    group = as_factor(c(
        rep("Peer Conformity", 2), "Market Competition")),
      conf.low.95 = estimate - std.error * 1.96,
      conf.high.95 = estimate + std.error * 1.96,
      .width = 0.95,
      conf.low.90 = estimate - std.error * 1.64,
      conf.high.90 = estimate + std.error * 1.64,
      .width = 0.9) %>%
    plot_ols_estimates(., "label") +
    my_theme() +
    ggtitle("")  +
    theme(legend.position = "right") +
    xlab("Changes in Predicted Probabilities") 

ggsave("plot/ukraine_ceo.pdf", width = 8, height = 8/1.618)


# -----------------------------------
# industry heterogeneity -- Figure D.6
# -----------------------------------

INDUSTRIES <- unique(u_sample$industry_group)

covariate_ind = c("years_emp", "position", "college", "income") 
covariate_firm = c("TKY_OSK", "estab_year", "foreign_owned", "capital")
pretreat_covs = c("sanction", "sanction_impact", "sanction_second")

industry_models = c()

for (k in c(1:length(INDUSTRIES))){

  df_industry <- u_sample %>% filter(industry_group != INDUSTRIES[[k]])
  
  collect_models_ind = c()

  for(i in 1:3){
    treatment = treatment_vec[i]
    control = control_vec[i]
    
    run_lm(
        data = df_industry,
        treatment = treatment,
        control = control,
        outcome = "u_stop",
        covs = c(pretreat_covs, covariate_firm, covariate_ind),
        covs_include = TRUE)  -> model_ind
    
    collect_models_ind[[i]] = model_ind
  }

  # get results
  bind_rows(
    collect_models_ind[[3]] %>%
        get_data("US withdrawal"),
    collect_models_ind[[2]] %>%
        get_data("Multiple withdrawal"),
    collect_models_ind[[1]] %>%
        get_data("China stays")) %>%
    mutate(
        industry = paste(INDUSTRIES[[k]], "(excluded)")) -> collect_industry
    
  industry_models[[k]] = collect_industry
  
}


est_industry_comb <- do.call(rbind, industry_models)


est_industry_comb %>% 
  mutate(
      label = factor(label, levels = c("China stays", "US withdrawal", "Multiple withdrawal")),
      group  = as_factor(industry),
      group = recode(group, "textile_furniture (excluded)" = "w/o Textile and Furniture", "food_beverage (excluded)" = "w/o Food and Beverage", "chemical_metal (excluded)" = "w/o Chemical and Metal", "metal (excluded)" = "w/o Metal", "machinery (excluded)" = "w/o Machinery", "construction_mining (excluded)" = "w/o Construction and Mining", "transportation (excluded)" = "w/o Transportation", "others (excluded)" = "w/o Others"),
      conf.low.95 = estimate - std.error * 1.96,
      conf.high.95 = estimate + std.error * 1.96,
      .width = 0.95,
      conf.low.90 = estimate - std.error * 1.64,
      conf.high.90 = estimate + std.error * 1.64,
      .width = 0.9) %>% 
    plot_ols_estimates(., "label") +
    my_theme() +
  xlab("Estimates") + 
  ggtitle("") +
  geom_vline(xintercept = 0, color = '#6f6f6f', linetype = "dashed") +
  ylab("") + 
  theme(legend.position = "right") 


ggsave("plot/ukraine_industry.pdf", width = 10, height = 10)


# ----------------------
# ordered logit -- Figure D.7
# ----------------------

covariate_ind = c("years_emp", "position", "college", "income")
covariate_firm = c("industry_group", "employee", "TKY_OSK", "estab_year", "foreign_owned")
pretreat_covs = c("sanction_impact")

# convert to factor
u_sample %>%
  mutate(u_stop_3level = as.factor(u_stop_3level)) -> u_sample

# China withdrawal
set.seed(02138)
run_simulation(
  8, 1500, bootstrap_orderlogit, 
  data = u_sample,
  treatment = "u_treatment_firm_china",
  control = "u_treatment_firm_multiple",
  outcome = "u_stop_3level",
  covariate_ind = covariate_ind,
  covariate_firm = covariate_firm,
  pretreat_covs = pretreat_covs,
  covs_include = TRUE) -> est1_china

# Multiple
run_simulation(
  8, 1500, bootstrap_orderlogit, 
  data = u_sample,
  treatment = "u_treatment_firm_multiple",
  control = "u_treatment_firm_us",
  outcome = "u_stop_3level",
  covariate_ind = covariate_ind,
  covariate_firm = covariate_firm,
  pretreat_covs = pretreat_covs,
  covs_include = TRUE) -> est1_multi


# US
run_simulation(
  8, 1500, bootstrap_orderlogit, 
  data = u_sample,
  treatment = "u_treatment_firm_us",
  control = "u_treatment_control",
  outcome = "u_stop_3level",
  covariate_ind = covariate_ind,
  covariate_firm = covariate_firm,
  pretreat_covs = pretreat_covs,
  covs_include = TRUE) -> est1_us

  # Figure 2
density_data_fig2 <- bind_rows(
  est1_china %>% get_filtered_data("China stays"),
  est1_multi %>% get_filtered_data("Multiple withdrawal"),
  est1_us %>% get_filtered_data("US withdrawal")
)

density_data_fig2 %>%
  filter(names == "level3") %>% 
  group_by(label) %>%
  mutate(
    conf.low.95 = quantile(estimate, 0.025),
    conf.high.95 = quantile(estimate, 0.975),
    conf.low.90 = quantile(estimate, 0.05),
    conf.high.90 = quantile(estimate, 0.95),
    est = mean(estimate),
    cov = "cov",
    label = factor(label, levels = c("China stays", "US withdrawal", "Multiple withdrawal"))) %>%
  dplyr::select(-estimate) %>% distinct() %>%
  plot_orderlogit_estimates(
    ., term = "label", estimate = "est") +
    geom_vline(xintercept = 0, color = '#6f6f6f', linetype = "dashed") +
    my_theme() + ylab("") + xlab("Change in Predicted Probabilities") +
    scale_x_continuous(limits = c(-0.2, 0.2)) +
    theme(legend.position= "none") +
    ggtitle("Effects on Support for Withdrawal from Russia") 

ggsave("plot/ukraine_orderLogit.pdf", width = 7, height = 5)


# ----------------------
# Japanese government policy and public opinion -- Figure D.8
# ----------------------


covariate_ind = c("years_emp","position","college","income")
covariate_firm = c("industry_group","employee","TKY_OSK","capital","sales")
pretreat_covs = c("sanction", "sanction_impact", "sanction_second")

# run models
run_logit_all(
  "info_japan_public", 
  "info_japan_public", 
  "info_japan_public") -> out_jp_public

# run models
run_logit_all(
  "info_japan_support", 
  "info_japan_support", 
  "info_japan_support") -> out_jp_gov


# US
bind_rows(
  bind_rows(out_jp_public$est1, out_jp_public$est2) %>% 
    filter(term == "u_treatment_firm_us") %>%
    mutate(model = recode_factor(model, "US withdrawal" = "Japanese public opinion")),
  bind_rows(out_jp_gov$est1, out_jp_gov$est2) %>% 
    filter(term == "u_treatment_firm_us") %>%
    mutate(model = recode_factor(model, "US withdrawal" = "Japense government policy"))
  ) %>%
  mutate(
    cov = factor(model, levels = c( "Japanese government policy", "Japanese public opinion"))) %>%
  plot_logit_estimates(., "model") +
  theme(legend.position= "none") +
  ggtitle("US withdrawal") +
  xlab("Change in Probabilities to\nSeek More Information") +
  theme(axis.title.y=element_blank(),
        axis.ticks.y=element_blank(),
        axis.text.y=element_blank(),
        plot.margin = margin(0.5, 0.5, 0.5, 0.8, "cm")) +
  scale_x_continuous(limits= c(-1, 1)) -> plot_jp_us

# Multiple 
bind_rows(
  bind_rows(out_jp_public$est1, out_jp_public$est2) %>% 
    filter(term == "u_treatment_firm_multiple") %>%
    mutate(model = recode_factor(model, "Multiple withdrawal" = "Japanese public opinion")),
  bind_rows(out_jp_gov$est1, out_jp_gov$est2) %>% 
    filter(term == "u_treatment_firm_multiple") %>%
    mutate(model = recode_factor(model, "Multiple withdrawal" = "Japense government policy"))
  ) %>%
  mutate(
    cov = factor(model, levels = c( "Japanese government policy", "Japanese public opinion"))) %>%
  plot_logit_estimates(., "model") +
  theme(legend.position= "none") +
  ggtitle("Multiple withdrawal") +
  xlab("") +
  theme(axis.title.y=element_blank(),
        axis.ticks.y=element_blank(),
        plot.margin = margin(0.5, 0.5, 0.5, 0.8, "cm")) +
  scale_x_continuous(limits= c(-1, 1)) -> plot_jp_multi

# China
bind_rows(
  bind_rows(out_jp_public$est1, out_jp_public$est2) %>% 
    filter(term == "u_treatment_firm_china") %>%
    mutate(model = recode_factor(model, "China stays" = "Japanese public opinion")),
  bind_rows(out_jp_gov$est1, out_jp_gov$est2) %>% 
    filter(term == "u_treatment_firm_china") %>%
    mutate(model = recode_factor(model, "China stays" = "Japense government policy"))
  ) %>%
  mutate(
    cov = factor(model, levels = c( "Japanese government policy", "Japanese public opinion"))) %>%
  plot_logit_estimates(., "model") +
  theme(legend.position= "none") +
  ggtitle("China stays") +
  xlab("") +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        plot.margin = margin(0.5, 0.5, 0.5, 0.8, "cm")) +
  scale_x_continuous(limits= c(-1, 1)) -> plot_jp_china

plot_jp_multi + plot_jp_us + plot_jp_china -> plot_jp

plot_jp +
  plot_annotation(
    title = 'Heterogeneous Effects on Seeking Information',
    theme = theme(plot.title = element_text(size = 20, hjust = 0.7, face="bold")))

ggsave("plot/ukraine_H7-4_jp_combined.pdf", width = 12, height = 6)

